/* peq.reajs
 * 
 * Copyright (c) 2012, Michael Gruhn
 * All rights reserved.
 * 
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted.
 * 
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY, FITNESS AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
 * PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
 * THIS SOFTWARE.
 */

/* Parametric EQ after [1].
 */
desc: Parametric EQ

slider1:12000<20,22000,0.001>Frequency (Hz)
slider2:0<-12,120,0.000001>Gain (dB)
slider3:1<0,10,0.0001>Bandwidth (Oct)
slider5:0<0,1,1>-bypass
slider9:5<0,7,{lowpass,highpass,bandpass,notch,allpass,peak,lowshelf,highshelf}>Selector
slider11:0<-240,24,0.0001>drive
slider13:1<1,200,1>-order
slider14:500<100,100000,1>plotsize
slider15:1<1,1000,1>bands
slider16:12<0.01,240,0.000001>gain range
slider17:0<0,1,{lookup table, actual magnitude}>-mode
//slider9:0<-0.000001,0.000001,0.000000000000001>-gain
@init
b0=0;
b1=1000;
b2=2000;
a0=3000;
a1=4000;
a2=5000;

bufferl=6600;
bufferr=6800;
buftemp=7000;

freqplot=100000;
phaseplot=300000;
function tanh(x)
     (
x = exp(2*x); 
(x - 1) / (x + 1);  //(ex - e-x)/(ex + e-x)
);
function atanh(x)
     (
       (1/2*log((1+x)/(1-x)));
     );
     function sinh(x)
     (
((e^x) - (e^-x))/2;
);

function cosh(x)
     (
((e^x) + (e^-x))/2;
);

function tanh(x,a)
     (
(($e^(x*a)) - ($e^(-x*(1-a))))/(($e^(x*a)) + ($e^(-x*(1-a))));  //(ex - e-x)/(ex + e-x)
);

function coth(x)
     (
((e^x) + (e^-x))/((e^x) - (e^-x));
);

function sech(x)
     (
2/((e^x) + (e^-x));
);

function csch(x)
     (
2/((e^x) - (e^-x))/2;
);
function bezier(y1,y2,y3,y4,t) (
        y1 + 3*t*(y2-y1) + 3*(t^2)*(y1+y3-(2*y2)) + (t^3)*(y4-y1+3*y2-3*y3);
);    
@slider
c = 8.65617025;
plotsize=floor(slider14);
drive = pow(10,slider11/20);
drive2 = pow(10,slider10/20);
driveb = pow(10,slider11/20);
drive2b = pow(10,slider12/20);
invdrive = pow(10,(-196)/20);
A = pow(10,(slider2/40));
q=1/slider3;

 w0 = (2 * $pi * (slider1)) / srate;
 w02 = (2 * $pi * (-slider1)) / srate;
        alpha = sin(w0) / (2 * q);//*1/sqrt(3);
        alpha2 = alpha;//*A+alpha*(1-a);
slider9 == 0 ? (
            b0[passes] = (1.0 - cos(w0)) / 2;
            b1[passes] = 1.0 - cos(w0);
            b2[passes] = (1.0 - cos(w0)) / 2;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2* cos(w0);
            a2[passes] = 1.0 - alpha;
        );
 slider9 == 1 ? (
            b0[passes] = (1.0 + cos(w0)) / 2;
            b1[passes] = -(1.0 + cos(w0));
            b2[passes] = (1.0 + cos(w0)) / 2;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha;
        );
    slider9 == 2 ? (
            b0[passes] = alpha;
            b1[passes] = 0.0;
            b2[passes] = -alpha;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha;
       );
slider9 == 3 ? (
            b0[passes] = 1.0;
            b1[passes] = -2 * cos(w0);
            b2[passes] = 1.0;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha;
        );
slider9 == 4 ? (
            b0[passes] = 1.0+ (alpha2);
            b1[passes] = -2 * cos(w0);
            b2[passes] = 1.0- (alpha2);
            a0[passes] = 1.0+ (alpha2);
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0- (alpha2);
       );
slider9 == 5 ? (
            b0[passes] = 1.0 + alpha * A;
            b1[passes] = -2 * cos(w0);
            b2[passes] = 1.0 - alpha * A;
            a0[passes] = 1.0 + alpha / A;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha / A;
            d0=1.0 + alpha * A;
            d1=-2*cos(w0);
            d2=1.0 - alpha * A;
       );
slider9 == 6 ? (
            b0[passes] = A * (((A + 1.0) - (A - 1.0) * cos(w0)) + 2 * sqrt(A) * alpha);
            b1[passes] = 2 * A * (A - 1.0 - (A + 1.0) * cos(w0));
            b2[passes] = A * ((A + 1.0) - (A - 1.0) * cos(w0) - 2 * sqrt(A) * alpha);
            a0[passes] = A + 1.0 + (A - 1.0) * cos(w0) + 2 * sqrt(A) * alpha;
            a1[passes] = -2 * ((A - 1.0) + (A + 1.0) * cos(w0));
            a2[passes] = (A + 1.0 + (A - 1.0) * cos(w0)) - 2 * sqrt(A) * alpha;
       );
 slider9 == 7 ? (
            b0[passes] = A * (A + 1.0 + (A - 1.0) * cos(w0) + 2 * sqrt(A) * alpha);
            b1[passes] = -2 * A * ((A - 1.0) + (A + 1.0) * cos(w0));
            b2[passes] = A * ((A + 1.0 + (A - 1.0) * cos(w0)) - 2 * sqrt(A) * alpha);
            a0[passes] = ((A + 1.0) - (A - 1.0) * cos(w0)) + 2 * sqrt(A) * alpha;
            a1[passes] = 2 * (A - 1.0 - (A + 1.0) * cos(w0));
            a2[passes] = (A + 1.0) - (A - 1.0) * cos(w0) - 2 * sqrt(A) * alpha;
        );
        /*
        scaler = srate/2;
        tmp=a0*scaler;
        a0=(1.0+tmp)/(1.0-tmp);
        tmp=a1*scaler;
        a1=(1.0+tmp)/(1.0-tmp);
        tmp=a2*scaler;
        a2=(1.0+tmp)/(1.0-tmp);
        tmp=b0*scaler;
        b0=(1.0+tmp)/(1.0-tmp);
        tmp=b1*scaler;
        b1=(1.0+tmp)/(1.0-tmp);
        tmp=b2*scaler;
        b2=(1.0+tmp)/(1.0-tmp);
       */
       /*
        c1[passes] = b0[passes] / a0[passes];
        c2[passes]= b1[passes] / a0[passes];
        c3[passes]= b2[passes] / a0[passes];
        c4[passes]= a1[passes]/ a0[passes];
        c5[passes]=a2[passes]/ a0[passes];
        c6[passes]= d0 / a0[passes];
        c7[passes]= d1/ a0[passes];
        c8[passes]=d2/ a0[passes];
        */
        //b0-=b0b;
       //     b1-=b1b;
       //    b2-=b2b;
       //    a0-=a0b;
      //     a1-=a1b;
       //     a2-=a2b;
  
redraw = 25;

@sample

spl0*=drive;
spl1*=drive;


slider17 == 0 ? (

bufferl[i]=spl0;
i2=0;
loop(32,
buftemp[i2]=bufferl[(i+i2)%32];
i2+=1;
);
memset(buftemp+32,0,32);
fft(buftemp,32);
fft_permute(buftemp,32);
mean=0;
oldphase=phasemean;
phasemean=0;
i2=0;
maxfreq=0;
loop(32,
//maxfreq=max(maxfreq,abs(buftemp[i2]/64));
mean+=sqrt(buftemp[i2]^2+buftemp[i2+fftsize]^2)*(i2/30);
//phasemean+=atan2(buftemp[i2+fftsize],buftemp[i2]);
i2+=1;
);
//mean/=abs(spl0);
//mean/=62;
mean/=32;
mean=sin(mean*2*$pi);
//phasemean/=32*$pi;
//phasemean=abs(phasemean);
//mean/=62;
//mean=(mean);
//mean/=maxfreq;
//phase=abs(atan(phasemean))/($pi);//abs(atan2(phasemean,oldphase))/(2*$pi);
//mean=abs(mean);
//mean/=abs(spl0);
//phasemean/=128;
//phasemean=abs(phasemean);
//test3=sin(mean*2*$pi);
//test2=log(sqrt(mean));
test=max(0,min(plotsize,(mean)*plotsize));
//maxtest=0;
maxtest=max(maxtest,test);
interpolate=test-floor(test);//abs((freqplot[max(0,(floor(test-1)))]-freqplot[min(plotsize,ceil(test+1))]))/sqrt(freqplot[max(0,(floor(test-1)))]^2+freqplot[min(plotsize,ceil(test+1))]^2);
freqtemp=bezier(freqplot[max(0,(floor(test-1)))],freqplot[floor(test)],freqplot[min(plotsize,ceil(test))],freqplot[min(plotsize,ceil(test+1))],interpolate);
phasetemp=bezier(phaseplot[max(0,(floor(test-1)))],phaseplot[floor(test)],phaseplot[min(plotsize,ceil(test))],phaseplot[min(plotsize,ceil(test+1))],interpolate);
phasetemp=abs(phasetemp)/($pi);
//freqtemp*=test2;
//gain=max(0.0001,min(100,freqtemp));
gain=exp(freqtemp*(phasetemp)/c);
gain=atan(gain)/($pi/2);
//gain/=10^(3/20);
//foo=(gain)/200*$pi;
//bar = sin((gain)/200*$pi);
spl0 =tanh(spl0,gain);//+spl0*(1-gain);i+=1;
//spl0=tanh(spl0,gain);//*(phase)+spl1*(1-phase);//*(1-phasemean)+(spl0*(phasemean));
bufferr[i]=spl1;
i2=0;
loop(32,
buftemp[i2]=bufferr[(i+i2)%32];
i2+=1;
);
memset(buftemp+32,0,32);
fft(buftemp,32);
fft_permute(buftemp,32);
mean=0;
oldphase=phasemean;
phasemean=0;
i2=0;
maxfreq=0;
loop(32,
//maxfreq=max(maxfreq,abs(buftemp[i2]/64));
mean+=sqrt(buftemp[i2]^2+buftemp[i2+fftsize]^2)*(i2/30);
//phasemean+=atan2(buftemp[i2+fftsize],buftemp[i2]);
i2+=1;
);
//mean/=abs(spl0);
//mean/=62;
mean/=32;
mean=sin(mean*2*$pi);
//phasemean/=32*$pi;
//phasemean=abs(phasemean);
//mean/=62;
//mean=(mean);
//mean/=maxfreq;
//phase=abs(atan(phasemean))/($pi);//abs(atan2(phasemean,oldphase))/(2*$pi);
//mean=abs(mean);
//mean/=abs(spl0);
//phasemean/=128;
//phasemean=abs(phasemean);
//test3=sin(mean*2*$pi);
//test2=log(sqrt(mean));
test=max(0,min(plotsize,(mean)*plotsize));
//maxtest=0;
maxtest=max(maxtest,test);
interpolate=test-floor(test);//abs((freqplot[max(0,(floor(test-1)))]-freqplot[min(plotsize,ceil(test+1))]))/sqrt(freqplot[max(0,(floor(test-1)))]^2+freqplot[min(plotsize,ceil(test+1))]^2);
freqtemp=bezier(freqplot[max(0,(floor(test-1)))],freqplot[floor(test)],freqplot[min(plotsize,ceil(test))],freqplot[min(plotsize,ceil(test+1))],interpolate);
phasetemp=bezier(phaseplot[max(0,(floor(test-1)))],phaseplot[floor(test)],phaseplot[min(plotsize,ceil(test))],phaseplot[min(plotsize,ceil(test+1))],interpolate);
phasetemp=abs(phasetemp)/($pi);
//phasetemp*=10^(36/20);
//freqtemp*=test2;
//gain=max(0.0001,min(100,freqtemp));
gain=exp(freqtemp*(phasetemp)/c);
gain=atan(gain)/($pi/2);
//gain/=10^(3/20);
//foo=(gain)/200*$pi;
//bar = sin((gain)/200*$pi);
spl1 =tanh(spl1,gain);//+spl0*(1-gain);i+=1;
//spl1*=2;
);
i+=1;
i >= 32 ? i=0;
//spl0/=drive;
//spl1/=drive;
/*
i=0;
i2=0;
loop(2,
x0[i]=spl(i%2);
//tempy0=c1*x0 + c2*x1 + c3*x2 - c4*y1 - c5*y2;
y0[i] = c1[i2]*x0[i] + c2[i2]*x1[i] + c3[i2]*x2[i] - c4[i2]*y1[i] - c5[i2]*y2[i];

//phase = atan(y0/x0);
//spl0 = y0 * abs(phase-alpha)/$pi + (tanh(y0*drive))*(1-abs(phase-alpha)/$pi);
y2[i] = y1[i]; y1[i] = y0[i]; x2[i] = x1[i]; x1[i] = x0[i];
spl(i%2)=y0[i];


//x0r[i]=spl1;
//y0r[i] = c1*x0r[i] + c2*x1r[i] + c3*x2r[i] - c4*y1r[i] - c5*y2r[i];
/*
y0r = c1*x0r + c2*x1r + c3*x2r - c4*y1r - c5*y2r;
phase-alpha > 0 ? (
y0r >= 0 ? (
y0r = y0r * abs(phase-alpha)/$pi + (tanh(y0r*drive)/drive)*(1-abs(phase-alpha)/$pi);
) : (
y0r = y0r * abs(phase-alpha)/$pi + (tanh(y0r*driveb)/driveb)*(1-abs(phase-alpha)/$pi);
);
) : (
y0r >= 0 ? (
y0r = y0r * abs(phase-alpha)/$pi + (tanh(y0r*drive2)/drive2)*(1-abs(phase-alpha)/$pi);
) : (
y0r = y0r * abs(phase-alpha)/$pi + (tanh(y0r*drive2b)/drive2b)*(1-abs(phase-alpha)/$pi);
);
);
*/
//y2r[i] = y1r[i]; y1r[i] = y0r[i]; x2r[i] = x1r[i]; x1r[i] = x0r[i];
//spl0=y0;
//s//pl1=y0r[i];
/*
i%2 == 1 ? i2+=1;
i+=1;
);
*/
@gfx
mouse_wheel > 0 ? (
passes+=1;
slider2=0;
slider1=16000;
slider3=4;
);
mouse_cap & 1 ? (
tmp = (mouse_x/gfx_w)*srate/2;
 // g_w = 2*$pi*tmp/srate;
slider1 = min(22000,max(16,tmp));
slider2 = ((-(mouse_y/gfx_h)+0.5)*((slider16/2)*4));
);
mouse_cap & 2 == 0 ? (
rightmouseclick=0;
oldslider3=slider3;
);
mouse_cap & 2 ? (
rightmouseclick == 0 ? (
initx2=mouse_x;
rightmouseclick=1;
);
//slider3 = 20+(mouse_x/gfx_w)*19980;
slider3 = max(0.00001,(oldslider3+(((mouse_x-initx2)/gfx_w))*(oldslider3*2)));
);
mouse_wheel > 0 || mouse_cap > 0 ? (
A = pow(10,(slider2/40));
q=1/slider3;

 w0 = (2 * $pi * (slider1)) / srate;
 w02 = (2 * $pi * (-slider1)) / srate;
        alpha = sin(w0) / (2 * q)*1/sqrt(3);
        alpha2 = alpha;//*A+alpha*(1-a);
slider9 == 0 ? (
            b0[passes] = (1.0 - cos(w0)) / 2;
            b1[passes] = 1.0 - cos(w0);
            b2[passes] = (1.0 - cos(w0)) / 2;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2* cos(w0);
            a2[passes] = 1.0 - alpha;
        );
 slider9 == 1 ? (
            b0[passes] = (1.0 + cos(w0)) / 2;
            b1[passes] = -(1.0 + cos(w0));
            b2[passes] = (1.0 + cos(w0)) / 2;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha;
        );
    slider9 == 2 ? (
            b0[passes] = alpha;
            b1[passes] = 0.0;
            b2[passes] = -alpha;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha;
       );
slider9 == 3 ? (
            b0[passes] = 1.0;
            b1[passes] = -2 * cos(w0);
            b2[passes] = 1.0;
            a0[passes] = 1.0 + alpha;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha;
        );
slider9 == 4 ? (
            b0[passes] = 1.0+ (alpha2);
            b1[passes] = -2 * cos(w0);
            b2[passes] = 1.0- (alpha2);
            a0[passes] = 1.0+ (alpha2);
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0- (alpha2);
       );
slider9 == 5 ? (
            b0[passes] = 1.0 + alpha * A;
            b1[passes] = -2 * cos(w0);
            b2[passes] = 1.0 - alpha * A;
            a0[passes] = 1.0 + alpha / A;
            a1[passes] = -2 * cos(w0);
            a2[passes] = 1.0 - alpha / A;
            d0=1.0 + alpha * A;
            d1=-2*cos(w0);
            d2=1.0 - alpha * A;
       );
slider9 == 6 ? (
            b0[passes] = A * (((A + 1.0) - (A - 1.0) * cos(w0)) + 2 * sqrt(A) * alpha);
            b1[passes] = 2 * A * (A - 1.0 - (A + 1.0) * cos(w0));
            b2[passes] = A * ((A + 1.0) - (A - 1.0) * cos(w0) - 2 * sqrt(A) * alpha);
            a0[passes] = A + 1.0 + (A - 1.0) * cos(w0) + 2 * sqrt(A) * alpha;
            a1[passes] = -2 * ((A - 1.0) + (A + 1.0) * cos(w0));
            a2[passes] = (A + 1.0 + (A - 1.0) * cos(w0)) - 2 * sqrt(A) * alpha;
       );
 slider9 == 7 ? (
            b0[passes] = A * (A + 1.0 + (A - 1.0) * cos(w0) + 2 * sqrt(A) * alpha);
            b1[passes] = -2 * A * ((A - 1.0) + (A + 1.0) * cos(w0));
            b2[passes] = A * ((A + 1.0 + (A - 1.0) * cos(w0)) - 2 * sqrt(A) * alpha);
            a0[passes] = ((A + 1.0) - (A - 1.0) * cos(w0)) + 2 * sqrt(A) * alpha;
            a1[passes] = 2 * (A - 1.0 - (A + 1.0) * cos(w0));
            a2[passes] = (A + 1.0) - (A - 1.0) * cos(w0) - 2 * sqrt(A) * alpha;
        );
        redraw=25;
);
mouse_wheel < 0 ? (
passes-=1;
);
mouse_wheel = 0;
(redraw) >= 25 ? (
redraw = 0;

gfx_r=gfx_g=gfx_b=0; gfx_a=1;
gfx_x=gfx_y=0;
gfx_rectto(gfx_w,gfx_h);
/*
gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_y = 0; gfx_x = (log(20)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(30)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(40)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(50)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(60)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(70)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(80)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(90)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_y = 0; gfx_x = (log(100)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_y = 0; gfx_x = (log(200)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(300)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(400)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(500)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(600)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(700)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(800)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(900)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_y = 0; gfx_x = (log(1000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_y = 0; gfx_x = (log(2000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(3000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(4000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(5000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(6000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(7000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(8000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
gfx_y = 0; gfx_x = (log(9000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_y = 0; gfx_x = (log(10000)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);

gfx_r=0.5; gfx_g=gfx_b=0; gfx_a=1;
gfx_y = 0; gfx_x = (log(srate/2)/log(10)-log(20)/log(10))/3*gfx_w; gfx_lineto(gfx_x, gfx_h, 0);
*/
gfx_r=gfx_g=gfx_b=1; gfx_a=1;
gfx_x = 0; gfx_y = gfx_h - (0 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
/*
gfx_r=gfx_g=gfx_b=0.5; gfx_a=1;
gfx_x = 0; gfx_y = gfx_h - (6 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (-6 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);

gfx_r=gfx_g=gfx_b=0.25; gfx_a=1;
gfx_x = 0; gfx_y = gfx_h - (3 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (-3 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (9 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
gfx_x = 0; gfx_y = gfx_h - (-9 + 12)/24 * gfx_h; gfx_lineto(gfx_w,gfx_y,0);
*/

g_x = 0;
gfx_x = 0;
gfx_y = gfx_h/2;
//plotsize=100;
testa=log(10);
testb=log(0);
i21=0;
loop(plotsize,
g_db=0;
i8=0;
ca=0;
sa=0;
cb=0;
sb=0;
loop(passes+1,
//  g_w = 2*$pi*(((g_x/gfx_w)^g_s*19980+20)/srate);
  tmp = (g_x/gfx_w)*srate;//10^(g_x/gfx_w*4);
  g_w = $pi*(tmp)/srate;
  g_phi = 4*sin(g_w/2)^2;
  g_db += ((10*log( (b0[i8]+b1[i8]+b2[i8])^2 + ( b0[i8]*b2[i8]*g_phi - (b1[i8]*(b0[i8]+b2[i8]) + 4*b0[i8]*b2[i8]) )*g_phi )/log(10) - 10*log( (a0[i8]+a1[i8]+a2[i8])^2 + ( a0[i8]*a2[i8]*g_phi - (a1[i8]*(a0[i8]+a2[i8]) + 4*a0[i8]*a2[i8]) )*g_phi )/log(10)));
  //g_db = (((g_db)));
  //g_db = g_db2
  //g_db = 10*log( (b0+b1+b2)^2 + ( b0*b2*g_phi - (b1*(b0+b2) + 4*b0*b2) )*g_phi )/log(10)
 // - 10*log( (a0+a1+a2)^2 + ( a0*a2*g_phi - (a1*(a0+a2) + 4*a0*a2) )*g_phi )/log(10);
 ca += 1 + a1[i8]*cos(g_w) + a2[i8]*cos(g_w);
 sa += a1[i8]*sin(g_w) + a2[i8]*sin(g_w);
 cb += b0[i8] + b1[i8]*cos(g_w) + b2[i8]*cos(g_w);
 sb += b1[i8]*sin(g_w) + b2[i8]*sin(g_w);
 
 i8+=1;
  );
 // g_db=min(100,max(0,g_db));
  g_y = gfx_h - (g_db + slider16*2)/(slider16*4) * gfx_h;
  //g_db=abs(g_db);
  gfx_r=0;gfx_g=1;gfx_b=0; gfx_a=1;
  
  gfx_lineto(g_x+1,g_y,0);
  
  gfx_x=g_x;
  gfx_y=g_y;
  phaseplot[i21] = atan2(sb,cb)-atan2(sa,ca);
  phaseplot[i21+1]=phaseplot[i21];
  freqplot[i21]=g_db;
  freqplot[i21+1]=g_db;
  i21+=2;
  g_x+=(1/(plotsize-1))*gfx_w;
  tempy=g_x;
);
//freqplot[i21]=g_db;
memset(freqplot+plotsize,g_db,plotsize);
gfx_r=0;gfx_g=0.5;gfx_b=1; gfx_a=1;
g_x = 0;
gfx_x = 0;
gfx_y = gfx_h/2;
);

/* 
 * References:
 * [1] Digital Parametric Equalizer Design With Prescribed Nyquist-Frequency Gain
 *     Sophocles J. Orfanidis
 */
slider15=passes;
@serialize
file_mem(0,0,10000);
file_var(0,passes);

